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Dynamical mean-field theory is generalized to solve the nonequilibrium Keldysh boundary prob- 
lem: a system is started in equilibrium at a temperature T = 0.1, a uniform electric field is turned 
on at t = 0, and the system is monitored as it approaches the steady state. The focus here is on 
the Bloch oscillations of the current and how they decay after their initial appearance near t — 0. 
The system is evolved out to the longest time allowed by our computational resources — in most 
cases we are unable to reach the steady state. The strongly correlated material is modeled by the 
spinless Falicov-Kimball model at half-filling on a hypercubic lattice in d = oo dimensions, which 
has a Mott-like metal-insulator transition at U = \/2- The computational algorithm employed is 
highly efficient, parallelizes well, and scales to thousands of processors. For strong fields, we find 
beats develop with a period of 2n/U, while for strong interactions, the Bloch oscillations are sharply 
damped and become quite irregular in time. 



PACS numbers: 71.27. 



71.10.Fd, 71.45.Gm, 72.20.Ht 



I 

O 

o 



> 

On 
O 



O 



X 



I. INTRODUCTION 



was 



mean-field theory (DMFT) 
in 1989 as a new technique to solve 



Dynamical 
introduced^ 

the many-body problem. Originally, DMFT was used 
to find the charge-density-wave transition temperature 
of the spinless Falicov-Kimball model- at half-filling. 
Since then, it has been employed to solve nearly all 
model solid-state Hamiltonians^ in equilibrium. It was 
recently generalized to treat nonequilibrium situations^. 
This contribution is a sequel to that work, where we 
examine the transient evolution of the Bloch oscillations 
of the current that appear in a system that starts in 
equilibrium and has a uniform electric field turned on at 
time t = — the so-called Keldysh boundary problem^. 

Bloch^ and Zener''' showed that noninteracting elec- 
trons (on a lattice) undergo an oscillatory motion when 
placed in a uniform static electric field (oriented along a 
symmetry direction of the Brillouin zone), because the 
electron wavevector, which evolves linearly in time under 
the driving force from the electric field, is Bragg reflected 
whenever it reaches a Brillouin zone boundary (since the 
wavevector will periodically trace out an identical orbit 
when the electric field is oriented in a symmetry direc- 
tion). But in conventional metals, Bloch oscillations have 
never been seen, because the electron relaxation time is 
so short, the electrons are scattered before they reach 
the zone boundary and Bragg reflect. Bloch oscillations 
have been observed in semiconducting heterostructures^, 
Josephson junctions^, and cold-atom systemsiS. In this 
contribution, we use nonequilibrium DMFT to solve the 
Keldysh boundary problem and analyze how the Bloch 
oscillations are quenched in a strongly correlated material 
as the scattering is increased. In particular, we describe 
the change in the behavior of the current as the system 
evolves from a metal to a Mott insulator with increasing 
U. 

The formalism for the nonequilibrium many-body 



problem was developed in the 1960s by Kadanoff and 
Baymii and Keldysh^. They determined the generic 
equations for the Green's functions and then solved them 
under different approximate conditions. Since then, the 
exact solution has remained one of the long-standing un- 
solved problems of many-body physics. Now, DMFT can 
be employed to solve this problem in the limit of large 
spatial dimensions. In this contribution, we describe, in 
detail, how to do this. This problem is quite complex, 
and necessarily, we need to limit our coverage to focus 
the discussion. We have chosen to concentrate on the 
phenomena of the quenching of Bloch oscillations, exam- 
ining how they change their character as one goes from 
a metal to a Mott insulator, and how they develop beats 
for large fields. 

The remainder of this contribution is organized as fol- 
lows: in Section II, we discuss the formalism, and derive, 
in detail, all of the relevant generalizations of DMFT 
to treat the nonequilibrium case as described by the 
Keldysh boundary problem. In Section III, we present 
our numerical results for extensive calculations of the cur- 
rent as a function of time on the hypercubic lattice. We 
will see a range of interesting behavior in both the high 
field limit and in the Mott insulator. We end with a 
discussion of the results and a conclusion in Section IV. 



II. FORMALISM 

The many-body formalism for nonequilibrium dy- 
namical mean-field theory is straightforward to develop 
within the Kadanoff-Baym-Keldysh approach-'--. Be- 
cause nonequilibrium problems are not time-translation 
invariant (indeed an external field is turned on at i = 0), 
we need to employ Green's functions that depend on 
two times. We work with the so-called contour-ordered 
Green's function, which is defined for any two time values 
that lie on the Kadanoff-Baym-Keldysh contour shown in 
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Fig. [TJ The system starts in equilibrium until time t = 
when a uniform electric field is turned on. The contour 
starts at some time before the field is turned on, runs 
out to a maximal time, then returns to the original time, 
and finally moves parallel to the negative imaginary axis 
a distance /3 (equal to the inverse of the temperature of 
the original equilibrium distribution) . The evolution for- 
ward, and then backward, in time is necessary because 
there is no a priori relationship between the quantum 
states at large positive times, and the original equilib- 
rium distribution, so we must evolve backward in time 
to reach a known equilibrium distribution. 
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FIG. 1: Kadanoff-Baym-Keldysh contour for the two-time 
Green's functions in the nonequilibrium case. We take the 
contour to run from tmin = —5 to tmax and back, and then 
extend downward parallel to the imaginary axis for a distance 
of p. The field is turned on at t = 0; i.e., the vector potential 
is nonzero only for positive times. 

We will be working with a spatially uniform electric 
field. We describe it via a spatially uniform vector po- 
tential (with a vanishing scalar potential) in the so-called 
Hamiltonian gauge (see below). In this case, we pre- 
serve the translational invariance of the system in posi- 
tion space, so we can describe results in either real space 
or momentum space. Working in the Heisenberg pic- 
ture, where all time dependence is in the operators, the 
contour-ordered Green's function (in real space) is de- 
fined by 

Cr^it^t') ^ -^TrT,e-'^«-c,(^)c](^')/2e„ (1) 

+i9,{t', i)Tre-'3«-c](t')c,(0/2e„ (2) 

where cj (Cj) are the electron creation (annihilation) op- 
erators for conduction electrons at site i, the subscript eq 
denotes the equilibrium Hamiltonian (chosen at any time 
prior to when the field is turned on) , and the equilibrium 
partition function satisfies Z^q = Tr exp[— /JTieq], with /3 
the inverse temperature. The symbol 9c{t,t') is the gen- 
eralization of the theta function to the contour, and it 
equals 1 if t lies after t' on the contour and it equals 
otherwise. The time-dependence of the fermionic oper- 
ators is expressed in the Heisenberg representation with 
the time-evolution taking place along the contour (the 
times t and t' are any two times on the three-branch 
contour). 

By choosing the time variables to lie on different 
pieces of the contour, we can determine different types 
of Green's functions^^. In particular, if both t and t' 



lie on the upper real branch, we have the time-ordered 
Green's function 

Gl^{t,t') - -zTr7ie-^^-c,(0c](t')/2e5, (3) 

if both t and t' are on the lower real branch, we have the 
anti-time-ordered Green's function 

Gl^(t,t') = -zTr7ie-^^-c,(t)ct(i')/2e„ (4) 

where the bar denotes anti-time-ordering. Similarly, 
when t lies on the upper real branch and t' lies on the 
lower real branch, we have the so-called lesser Green's 
function 



(5) 



and when t lies on the lower real branch and t' lies on the 
upper real branch, we have the so-called greater Green's 
function 



G>{t,t') 



-zTre-^«-q(t)4(t')/2e 



(6) 



Finally, when both times lie on the imaginary branch, we 
have the thermal Green's function (because the Hamilto- 
nian is the equilibrium Hamiltonian here). There also are 
mixed Green's functions, where one time argument lies 
on a real branch and the other on the imaginary branch. 
These Green's functions are required to properly deter- 
mine the transient response at short times. They also 
enter when we map the lattice problem onto the impu- 
rity problem in an additional time-dependent field. Note 
that we can similarly define corresponding momentum- 
dependent Green's functions in analogy to the real-space 
Green's functions defined above. 

In addition to the Green's functions that can be di- 
rectly extracted from the contour-ordered Green's func- 
tion, there are two other important Green's functions 
that can be defined. The retarded Green's function is 
defined for real times, and satisfies 

G^^{tX)=-^0{t-t')TTe-^'''^^{cM,c]{t')}+/Z,q, (7) 

and the advanced Green's function satisfies 

Gt^it,t') ^zeit' ~t)TTe-P^^^{cS),4it')U/Z,q, (8) 

where the curly brackets denote the anticommutator. 

In nonequilibrium physics problems, we focus primar- 
ily on two of these different Green's functions — the re- 
tarded Green's function, which determines the quantum- 
mechanical states of the system, and the lesser Green's 
function, which determines how the electrons are dis- 
tributed amongst those quantum states. It turns out 
that all of the other Green's functions defined with real 
time arguments can be constructed from these two "fun- 
damental" Green's functions. 

We will consider a uniform electric field turned on 
at time t = 0. We ignore all magnetic field and rela- 
tivistic effects and assume the field is always uniform in 
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space. Then we can describe the field via a spatially uni- 
form time-dependent vector potential in the Hamiltonian 
gauge, where the scalar potential vanishes 



(9) 



here we have set the speed of light equal to one. The 
noninteracting contour-ordered Green's function in mo- 
mentum space is 



Gl 



\t,t') = -iTrr,e-^«™CkW4(t')/2e,. 



q,non 7 



(10) 

where we use creation and annihilation operators in mo- 
mentum space, and evaluate all averages with respect to 
the equilibrium noninteracting Hamiltonian 



n. 



eq,non — ^ ^^)^\^^^^^ 



(11) 



with the bandstructure satisfying ek — 
linitj^oo cos(ki)/\/d. The nearest neigh- 

bor hopping is chosen to scaled as i = t* and we 
use t* as our energy unit (the lattice constant a is also 
set equal to 1). The electric field is introduced via the 
vector potential, which is chosen to be 



A{t) = -9{t)Et, 



(12) 



for a uniform field turned on at i — 0. The time- 
dependent Hamiltonian is then constructed via the 
Peierls' substitution^^, where k — > k — eA{t) in the 
bandstructure. Hence the time-dependent noninteract- 
ing Hamiltonian becomes 



(13) 



Note that the noninteracting Hamiltonian commutes 
with itself at all times, even though it is time-dependent. 
This result makes it easy to exactly solve for the nonin- 
teracting Green's functions in a fielc^ii^. First note that 
the time-dependent creation and annihilation operators 
satisfy 



(14) 



Ck(*) = exp l^-i J [ek+0(t)eEt - H j "^k; (15) 

which follows by directly integrating their equations of 
motion (the integration in time runs over the contour 
from the starting point at tmm to the desired time t be- 
cause we need the time arguments for any location on the 
contour; tmin — —5 for the calculations presented here). 
Substituting into the definition for the noninteracting 
contour-ordered Green's function in Eq. pop then pro- 
duces the noninteracting Green's function directly (the 



integral over time lies on the contour between the points 
t' and t) 



G^'"°"(t, t') = z[/(ek - m) - e,{t, t')]e'^(*-*') 



X exp 



(F)eEt 



(16) 



The Fermi-Dirac distribution /(ek — A^) — 1/[1 -I- 
exp(/3{ek — A*})] enters from the initial equilibrium dis- 
tribution of the electrons (prior to the field being turned 
on). 

The noninteracting Green's function is quite compli- 
cated for an electric field pointing in an arbitrary direc- 
tion. It turns out that the simplest problem to consider 
is one where the field points along a diagonal direction 
E = £^(1, 1, 1, 1). Then the noninteracting Green's 
function depends on two band energies: ek (already de- 
fined above) and ek = lim^^oo ^t* X^iLi sinki/\/d in the 
following way: 

G:;r"(t,i') = *[/(e-At)-0c(t,O]e^^(*-*') 

X exp(~ij^ di{[0{-t) + 9{t)cos{eEt)]e 



-0{t)esm{eEt)} 



(17) 



The two band energies are distributed with a joint density 
of states that is a Gaussian in each variablei^ 



p{e,e) = -exp(-e^ -e^). 

TT 



(18) 



It will turn out that all momentum-dependent quantities 
we are interested in will depend only on e and e, so we 
reduce the infinite-dimensional Brillouin zone to a two- 
dimensional energy space. Note that in equilibrium, all e 
dependence drops out, and we have all objects of interest 
depending only on the bandstructure e. 

The contour-ordered Green's function for the interact- 
ing many-body problem satisfies Dyson's equation (with 
all time integrals running over the contour) 

GiM = G^'"°"(i,t') (19) 

dtj^ dt'G^^'^°'^{t, t)Y.^(i, i')GUt', t'), 

which can be viewed as defining the contour-ordered self- 
energy. We have written out the Dyson equation explic- 
itly in momentum space, and we used the fact that the 
self-energy is local (and hence independent of momen- 
tum) in DMFT, a fact that is proven next. 

The argument that the self-energy is local in d ^ cxd is 
based on two facts — the expansion of the self-energy in 
powers of the hoppingii which show that the self-energy 
is local (in equilibrium) and the Langreth rulea^S, which 
show how to relate perturbation theory in equilibrium 
to perturbation theory for the nonequilibrium case, and 
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thereby shows that the nonequiUbrium self-energy is also 
local. This follows because the many-body perturba- 
tion theory diagrams are identical in structure for both 
the equilibrium and nonequilibrium perturbation theo- 
ries, so the power-counting analysis (in inverse powers 
of the spatial dimension d) of Metzner guarantees that 
the nonequilibrium self-energy remains local in DMFT. 
In other words, the nonequilibrium DMFT problem can 
be mapped onto an impurity problem in time-dependent 
fields, just like the equilibrium problem, except now the 
fields (the Green's functions and the self-energy) have two 
time arguments. 

The DMFT algorithm employed to solve for the tran- 
sient response of nonequilibrium systems is closely re- 
lated to the iterative approach used in equilibrium^^. We 
develop each of the steps of this algorithm next, and then 
describe the differences between the nonequilibrium cal- 
culations and the equilibrium ones. 

The first step is to generalize the Hilbert transform 
employed in equilibrium DMFT to the nonequilibrium 
case. Starting from the local self-energy 'E'^{t,t'), which 
is a continuous matrix operator with both time vari- 
ables running over the contour, we need to construct 
the local contour-ordered Green's function Gf^^{t,t') = 
^\^G^{t,t'). This is done by first expressing the mo- 
mentum dependence in terms of the two band energies 
e and e, and employing the Dyson equation in Eq. 
Hence, 



Gl,{t,t') ^ de dep{e,e) 



[i~g: 



E") 'G:;r" (i,t'), (20) 



where I denotes the identity matrix operator [equal to 
Sc{t,t') — the Dirac delta function on the contour, which 
satisfies J^dt'6c{t,t')f{t') = f{t)], is given in 

Eq. (|17p , and there are a number of implicit matrix mul- 
tiplications over the time variables (which range over the 
contour) on the right hand side. The integrand for each e 
and e point is found by performing two continuous matrix 
operator multiplications and one continuous matrix op- 
erator inversion. Hence, the Hilbert transform for equi- 
librium problems generalizes to a two-dimensional inte- 
gral of a matrix valued function that requires two matrix 
multiplies and one matrix inversion to determine the in- 
tegrand. 

Once the local Green's function has been determined, 
the next step is to extract the dynamical mean-field 
X'^{t,t'). This is determined by first finding the effective 
medium by using Dyson's equation and then extracting 
the dynamical mean field from the effective medium. In 
particular, we have the effective medium satisfying 



G^,{tX)=[{GW 



{t,t'), 



(21) 



and then the dynamical mean field becomes 

X^{t,t') = (tdt + ^,)S,{t,t')~(G^,)-\t,t') (22) 
= {id^ + fi)S,{t, t') - (GlX\t. t') + t')- 



Here ^^ is the derivative operator along the contour, so 
it equals +d/dt along the upper branch, —d/dt on the 
lower branch, and —id/dr along the negative imaginary 
time axis. Note that we calculate A"^ from [Gq]~^, so we 
do not normally compute Gq to find the dynamical mean 
field. 

Once the dynamical mean field has been found, we 
need to determine the impurity Green's function for the 
impurity problem that evolves in the presence of the dy- 
namical mean field (which is nonzero on the entire con- 
tour). For the Falicov-Kimball model one can immedi- 
ately write down the solution to this problem (at the mo- 
ment, we do not have any numerical techniques to solve 
this problem for other models in the time representation). 
The spinless Falicov-Kimball model- involves single-band 
conduction electrons hopping on a lattice, and localized 
electrons which do not move but do interact with the 
conduction electrons when they are in the same unit cell 
via a screened Coulomb interaction U . The equilibrium 
lattice Hamiltonian (in the absence of a field) is then 



Tie 



^vq! to) 



(23) 

Here, we have c\ (Cj) create (annihilate) a spinless con- 
duction electron at site i, Wi — Q ov 1 is the localized 
electron number operator at site i, and /i is the conduc- 
tion electron chemical potential. Although the Falicov- 
Kimball model is one of the simplest many-body physics 
models, it does have a Mott metal-insulator transition, so 
it is interesting to use this model to examine how Bloch 
oscillations are quenched due to strong electron corre- 
lations (but the model does not include any Zener tun- 
neling effects because there are no higher energy bands). 
The solution to the impurity problem can be found by 
solving the equations of motion for the contour-ordered 
Green's function (the procedure is essentially the same 
as in equilibrium except the Green's functions depend 
on two times now and the times run over the contour) 
resulting in 

G^,np(t,t') = {l-{w{))[{l^^^+^Ji)5,{t,t')-\{t,t')]-^ 

+ {wi) [{td^ + M - U)6c{t, t') - X{t, t')r' , 

(24) 

with {wi) = J^ii'^i)/-^ average localized electron 
filling. 

The Dyson equation in Eq. (pij) is then employed to ex- 
tract the impurity self-energy from the impurity Green's 
function and the effective medium. The algorithm is then 
iterated until it converges (we usually require the Green's 
functions to converge pointwise to better than one part 
in 10^ in order to end the calculation). 

In summary, the basic structure of the iterative ap- 
proach to solving the DMFT equations^ continues to 
hold. We start with a guess for the self-energy (which is 
usually chosen to be equal to the equilibrium self-energy) , 
then we sum the momentum-dependent Green's function 
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over the Brillouin zone to produce the local Green's func- 
tion. Next the dynamical mean-field for the impurity 
problem is extracted by using Dyson's equation for the 
local Green's function and self-energy, the impurity prob- 
lem is solved in the dynamical mean-field to produce the 
impurity Green's function, and Dyson's equation is used 
again to extract the impurity self-energy. In the self- 
consistent solution of the DMFT equations, the impurity 
self-energy will be equal to the lattice self-energy. If they 
are different, then the new lattice self-energy is taken to 
be equal to the new impurity self-energy, and the loop 
is iterated until it converges. The noncquilibrium algo- 
rithm includes the following modifications from the equi- 
librium algorithm: (i) the summation over the Brillouin 
zone now requires at least a double integral over two band 
energies; (ii) the Green's functions are described by con- 
tinuous matrix operators with time indices that run over 
the contour; and (iii) the impurity problem solver must 
be generalized to the noncquilibrium case. 

The objects we work with on the contour are contin- 
uous matrix operators. Unfortunately, there is no way 
to directly calculate with such objects on a computer. 
Instead, we need to first discretize the contour, and de- 
fine the continuous matrix operators as the limit where 
the discretization size goes to zero of discrete matrices, 
whose matrix elements are identified with the average 
value of the matrix operators within the corresponding 
discrete intervals on the contour. By performing calcula- 
tions with particular discretizations, and then taking the 
limit where the discretization goes to zero (via an ex- 
trapolation procedure), we can approximate the results 
for the continuous matrix operators. The contour is dis- 
cretized in the following manner: we choose a real-time 
spacing At which varies from 0.1 to 0.014, and we fix 
the spacing along the imaginary axis to At = O.li so 
our largest matrix is 5700 x 5700 and we evaluate inte- 
grals over the contour by discrete summations using the 
leftpoint rectangular integration rule. The matrix oper- 
ators are general complex matrices, which are manipu- 
lated using standard linear algebra packages (LAPACK 
and BLAS). 

More precisely, the discretization process involves Nt 
points on the upper real branch (ranging from t„iin to 
tmax — At), Nt points on the lower real branch (ranging 
from tmax to tmin + At), and 100 points along the imag- 
inary axis (ranging from tmin to tmin — i(3 + O.lz, with 
P = 10); hence At = {tmax — train) /Nt- The discrete time 
values on the contour are then 

tj - -t„,„ + (j - l)At, l<j<Nt, (25) 

= tmax - [j ~Nt- l)At, A^t + 1 < J < 2Nt, 
= tm,n - 0.1z(j - 2Nt - 1), 2Nt + l<j < 2Nt + 100, 

where we used the fact that the discretization along the 
imaginary axis is always fixed at At = 0.1 in our calcu- 
lations. To calculate integrals over the contour, we use 
a leftpoint rectangular integration rule for discretizing 



integrals. 

„ 2Ar, + 100 

/ dtfit) = ^»/(*»)' (26) 

where the weights satisfy 

Wj = At, l<j< Nt, 

= -At, Nt + l<j<2Nt, 

= -O.li, 2Nt + l<j < 2Nt + 100. (27) 

The leftpoint integration rule evaluates the function at 
the "earliest" point in the time interval that has been 
discretized for the quadrature rule (which is the left hand 
side of the interval when we are on the upper real branch; 
earliest is meant with regards to the sense that the con- 
tour is traversed). 

Note that the contour-ordered Green's function satis- 
fies a boundary condition where we identify the points 
tmin with tmin — One Can show from the definition 
of the contour-ordered Green's function in Eq. ([1]), and 
the invariance of the trace with respect to the ordering 
of operators that G^i{tmin,t') = ~G%{tmin - iP,t') and 
G%{t,tmin) = —G1^{t,tmin — This is similar to the 
antiperiodicity property that the thermal Green's func- 
tions satisfy, and the proof is identical. 

The delta function changes sign along the negative 
real time branch, and is imaginary along the last branch 
of the contour in order to satisfy the property that 
J^dt'Scit,t')f{t') = f{t). In addition, we find that the 
numerics work better if the definition of the delta func- 
tion is done via "point splitting" (when we calculate the 
inverse of a Green's function) so that the delta function 
does not lie on the diagonal, but rather on the first sub- 
diagonal matrix (in the limit as At — > it becomes a 
diagonal operator). Because we identify the times tmin 
and t„iin — if3, the point splitting approach to the defi- 
nition of the delta function allows us to incorporate the 
correct boundary condition into the definition of the dis- 
cretized delta function. Hence, we define the discretized 
delta function in terms of the quadrature weights, in the 
following way 

Sc{ti,tj) = ly^'^^J+i' for integration over j, (28) 

— — — for integration over j, (29) 
Wi-i 

where ti and tj are two points on the discretized con- 
tour as described in Eq. (|25p. and Wi are the quadra- 
ture weights described in Eq. ((27)) . We have a different 
formula for integration over the first variable versus in- 
tegration over the second variable because we are using 
the leftpoint quadrature rule. Note that the formulas in 
Eqs. dill) and ((29]) hold only when i ^ 1. When i = 1, 
the only nonzero matrix element for the discretized delta 
function is the l,j = 2Nt + 100 matrix element, and it 
has a sign change due to the boundary condition that 
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the Green's function satisfies. The discretization of the 
derivative of the delta function on the contour is even 
more comphcated. It is needed to determine the inverse 
of the effective medium operator for the impurity. The 
derivative is calculated by a two-point discretization that 
involves the diagonal and the first subdiagonal. Since all 
we need is the discrete representation of the operator 



[idf + fj]Sc{t, t'), we summarize the discretization of that 
operator as follows 

[idt + ^ASc{tJ,tk) = i^MjkT^, (30) 



with the matrix Mjk satisfying 



/I 

-l-iM^i 1 
-1 - iMpL 1 



M. 



jk 



1 + lAtfi \ 




~l + iAtn 1 

-^ + iAt^J. 1 



4- At^ 1 



(31) 



here At = 0.1. The top third of the matrix corre- 
sponds to the upper real branch, the middle third to the 
lower real branch and the bottom third to the imaginary 
branch. Note that the operator [id^ + ^]5c is the in- 
verse operator of the Green's function of a spinless elec- 
tron with a chemical potential fj,. Hence the determi- 
nant of this operator must equal the partition function 
of a spinless electron in a chemical potential /i, namely 
1 + exp[/3/i]. Taking the determinant of the matrix Mjk 
gives 

detM = l + {-lf^'+^-{l + ^At^l){-l-iAtnf'~^ 
X {-l + iAtfi)^'{-l- ATfi)^^ 
« l + {l + AT^i)^-+0{At^), (32) 

which becomes l-|-exp[/3/i] in the limit where At, At ~* 
{Nr is the number of discretization points on the imag- 
inary axis). This provides a check on our algebra, and 
shows the importance of the upper right hand matrix 
element of the operator, which is required to produce 
the correct boundary condition for the Green's function. 
This is also the reason why we chose to point-split the 
delta function when we defined it's discretized matrix 
operator. 

We also have to show how to discretize the continu- 
ous matrix operator multiplication and how to find the 
discretized approximation to the continuous matrix op- 
erator inverse. Matrix multiplication is discretized as fol- 
lows 

f dtA{t,t)B{t,t') ^Y,A{U,tk)WkB{tk,tj). (33) 



So we must multiply the columns (or the rows) of the 
discrete matrix by the corresponding quadrature weight 
factors. This can be done either to the matrix on the left 
(columns) or to the matrix on the right (rows). To cal- 
culate the inverse, we recall the definition of the inverse 
for the continuous matrix operator 

J^dtA{t,t)A-\t,t') = Sc{t,t'), (34) 

which discretizes to 

MU, tk)WkA-\tk, tj) = ^5^j, (35) 

k * 

where here we do not need to point-split the delta func- 
tion. Hence, the inverse of the matrix is found by invert- 
ing the matrix WiA{ti, tj)Wj, or, in other words, we must 
multiply the rows and the columns by the quadrature 
weights before using conventional linear algebra inversion 
routines to find the discretized version of the continuous 
matrix operator inverse. 

Having resolved the technical details of the discretiza- 
tion, which allows us to evaluate approximations to the 
continuous matrix operators on a computer, we now move 
to the last technical aspect, which is how we perform the 
numerical quadrature of a matrix-valued integral. Since 
the bandstructure energies e and e are both distributed 
with Gaussian weight functions, we employ Gaussian 
quadrature in each dimension to perform the integration. 
We found, by benchmarking the equilibrium solutior^, 
that averaging the results of two Gaussian quadratures 
with iV and TV -I- 1 points works better than choosing 
2N + 1 points for the quadrature. For most cases we 
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discuss here, we use iV = 54 and A'' + 1 = 55 points for 
the Gaussian quadrature routines, although some calcu- 
lations were performed with the more accurate iV = 100 
and + 1 = 101 routines. In the former case, we have a 
total of 5,941 quadrature points, while in the latter case 
it is 20,201 quadrature points. Since the calculation of 
each matrix in the integrand of the integral is indepen- 
dent of every other quadrature point, this part of the 
code is easily parallelized. The numerical quadrature re- 
quires two matrix multiplies and one matrix inversion for 
each quadrature point, and the integrand is a matrix. 

The calculation of the local Green's function from the 
local self-energy is the most computationally intensive 
part of the algorithm. Fortunately that part parallelizes 
well in the master-slave format, since we need to send the 
slave nodes the self-energy matrix, and then send them 
the energies e and e for the particular quadrature point, 
and the slaves calculate the integrand directly and accu- 
mulate the results locally. Then, once all bandenergies 
have been calculated, we use a recursive binary gather 
operational to accumulate the results and send them to 
the master efficiently. This works by dividing the slave 
nodes in half and having one half send their accumu- 
lated results to the other half, and having those slave 
nodes accumulate the results they received with the ones 
they had. A slave node that has sent its results to an- 
other slave node then becomes inactive. The procedure 
is repeated until only one slave node has all of the accu- 
mulated results, which is then sent to the master. The 
impurity solver which inputs the local Green's function 
for the lattice, and outputs the impurity self-energy, is a 
serial code, that cannot be parallelized, because the ma- 
trix operations need to all be performed in turn. These 
calculations are performed on the master node and they 
require a number of matrix inversions and matrix mul- 
tiplies, which are carried out with LAPACK and BLAS 
routines. We typically require between ten and one hun- 
dred iterations to reach convergence of the results (the 
total computer time for the calculations presented here 
was about 600,000 cpu-hours on a Cray XT3, 900,000 
cpu-hours on a SGI ALTIX, and 800,000 cpu-hours on a 
SUN OPTERON). Overall, the code is quite efficient. It 
has achieved over 65% of the peak speed on 2032 cores 
of an SGI ALTIX machine, and scales with near linear 
scaling up to many hundreds to a few thousand cores 
(the main issue that affects the linear scaling is the fact 
that part of the algorithm is serial in nature and does 
not scale). 

Once the Green's functions have converged, we calcu- 
late the current (in the Hamiltonian gauge) by evaluating 
the operator average 

(j(t)) =-ez^v[k + 0(i)eEt]G<(i,t). (36) 

k 

The velocity component is Vi (k) = r sin(ki)/2%/d, and 
all components of the current are equal when the field lies 
along the diagonal (hence we can replace by — e and 
change the sum over momentum to a two-dimensional 



integral over e and e). We also calculate the equal time 
retarded and lesser Green's functions and their first two 
derivatives (at equal time) and compare those results to 
the exact valued. In general, these "moments" are quite 
accurate as the step size is made smaller. 

We find, nevertheless, that the results for the current 
and for the moments usually need to be extrapolated to 
the limit At —^ 0. To do this, we use a Lagrange inter- 
polation formula to extrapolate the results to At = 0. 
Typically it is a quadratic extrapolation, requiring re- 
sults at three different values of At, but sometimes we 
use higher order extrapolants. This extrapolation proce- 
dure allows us to achieve quite accurate results for the 
moments, when compared to the exact values, and of the 
current (which must vanish when t < 0). Further de- 
tails of these numerical issues and of the accuracies are 
presented elsewhere2ii2^. 

We end our formalism discussion by defining the 
nonequilibrium many-body density of states. First we 
convert from the time variables t and t' to Wigner's aver- 
age T = {t + t')/2 and relative trei = t — t' time variables. 
Then the DOS is defined from the Fourier transform of 
the retarded Green's function with respect to the relative 
time. In other words, 

G«(T,a;) = J dtreie'^'-' (t + , T ^ ^-f^ 

(37) 

is the retarded Green's function at each average time, 
and the many-body density of states satisfies 

p°^^\T,to) = --lmG^iT,Lo). (38) 

TT 

We will use these relations to help analyze some of our 
results for the current. Note that the nonequilibrium 
density of states can actually become negative for finite 
average times — it is nonnegative before the field is turned 
on (in equilibrium) and also in the limit T — s- oo (the 
steady state). 

III. NUMERICAL RESULTS 



We produce numerical calculations of the nonequilib- 
rium current as a function of time for the case of half- 
filling, where the conduction electron and the localized 
electron fillings are each equal to 0.5. In equilibrium, this 
system has a metal-insulator transition at U = \pl. In 
the case where there is no scattering (C/ = 0), the Bloch 
oscillations continue forever. In the presence of scatter- 
ing (for metals), the Bloch oscillations maintain the same 
approximate "periodicity" , but the amplitude decays. As 
the interactions increase further, the oscillations become 
irregular. 

In Fig. [21 we plot the current for the noninteracting 
case, the case of a strongly scattering metal (J7 = 0.5, 
red), and the case of an anomalous metal (f/ = 1, green; 
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FIG. 2: Scaled nonequilibrium current for different values of 
U with E = 0.125. All of these cases are metals in equilib- 
rium. We used a quartic extrapolation formula with five At 
values for U = 0.5 (At = 0.1, 0.067, 0.05, 0.04, and 0.033), 
while we used a quadratic extrapolation formula with three 
At values for ?7 = 1 (At = 0.025, 0.02, and 0.0167). (Color 
on-line.) 



FIG. 3: Scaled nonequilibrium current for different values 
of U (0, 0.5, 1, and 1.5) with E = 0.25. All of these cases 
are metals in equilibrium except for [/ = 1.5 which is a near- 
critical Mott insulator. We used a quadratic extrapolation 
formula with three At values for U = 0.5 (At = 0.1, 0.067 
and 0.05, 0.04, and 0.033), (7 = 1 (At = 0.05, 0.04, and 
0.033), and U = 1.5 (At = 0.033, 0.025, and 0.02). (Color 
on-line.) 



anomalous in the sense that there is a dip in the many- 
body density of states near the chemical potential) for 
E = 0.125. The Bloch period in this case is IGtt « 50. 
The time range we can extend the calculations out to is 
tmax = 35, so we do not see one full Bloch oscillation in 
the time window. The initial temperature of the system 
satisfied /3 = 10, and the field is turned on at i = 0. 
We checked the errors of the extrapolated data against 
the known moment sum rules. The errors are largest at 
small times (in equilibrium, before the field is turned on) 
and improve for larger times. As a benchmark, we record 
the maximal error in the moments for times greater than 
t = 5. The exact moments are equal to 1 for the zeroth 
moment and 0.5 -f for the second moment. The 

U — 0.5 case has errors in the zeroth moment less than 
1% and errors of 2% for the second moment. The U — 1 
case has errors less than 1% for the zeroth moment and 
errors of 3% for the second moment. We find generically, 
that the calculations require smaller At values for smaller 
electric fields. 

In the strongly scattering metal, the Bloch oscillations 
appear to be simply damped, while in the anomalous 
metal, we start to see additional oscillations develop at 
short time, which appear to disappear at longer times. 
While we have clearly not reached the steady state yet, 
it does appear that the current is approaching a constant 
nonzero value for the anomalous metal at long times, as 
we expect for a driven correlated material. 

In Fig. [21 we plot the scaled current for E = 0.25 and 
four U values ranging from a ballistic metal to a near 
critical Mott insulator. The Bloch period is now approx- 
imately 25 units, so we fit about one and one half os- 
cillations in our time window. The behavior for these 
cases is similar in many respects to what we saw for 
the smaller field. The oscillations are damped and as 
the scattering strength increases we see some additional 
shorter period wiggles develop, especially for short times. 
The larger U values appear like they are approaching a 



steady-state limit (with a smaller value of the current 
than for E — 0.125, as expected for large fields^^), but 
it is also clear that the larger field strength increases the 
transient response and makes it more difficult to reach 
the steady state within our time window. We also find 
that our accuracy improves and we do not need to em- 
ploy as small values of At; this is a general trend that 
continues as we increase the electric field. The scaled mo- 
ments show the following accuracies: U = 0.5, less than 
1% error for the zeroth and second moments; C/ = 1, less 
than 1% error for the zeroth and second moments; and 
U — 1.5, less than 2% error for the zeroth moment and 
less than 4% error for the second moment. 

In Fig. 31 we plot the E = 0.5 current for normal met- 
als in panel (a) and for the anomalous metal and two 
Mott insulators in panel (b). Here we have nearly three 
full Bloch oscillations within our time window, and the 
metallic phases are behaving pretty much as we would 
expect — they show Bloch oscillations with the same pe- 
riod and are increasingly damped as U increases. Unfor- 
tunately, the larger field drives the system for a longer 
period of time, so we are unable to see the steady state 
limit emerging for any of these metallic cases. While 
we expect the steady state current to be a constant dc 
response, we cannot rule out the possibility of a small 
amplitude oscillatory response as well. The anomalous 
metal and insulating cases in panel (b) also display in- 
teresting behavior. For example, we see one prominent 
oscillation for the first half Bloch period and then the re- 
sponse is sharply damped afterwards, but the oscillations 
have no regularity to them, and as the system becomes 
more insulating we see additional shorter period oscilla- 
tions develop, similar to what occurred for smaller fields, 
but here they are more prominent. For longer times, the 
response is not yet reaching a steady state, although it 
is flattening out much faster. Notice that in some cases, 
the long-time current response is larger for more strongly 
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FIG. 4: Scaled nonequilibrium current for different values of 
U with E = 0.5. In panel (a), we show cases that are metals 
in equilibrium. In panel (b), we show one anomalous metal 
(U = 1) and two Mott insulators {U = 1.5 and U = 2); note 
the reduced scale for the vertical axis versus that in panel 
(a). We used a quadratic extrapolation formula for U < 1 
{At = 0.1, 0.067, and 0.05 for U = 0.125, U = 0.25 and 
U = 0.5; At = 0.067, 0.05, and 0.04 for f/ = 1), a cubic 
extrapolation formula for U = 1.5 {At = 0.067, 0.05, 0.04, 
and 0.033), and a quartic extrapolation formula for U — 2 
{At = 0.05, 0.04, 0.033, 0.025, and 0.02). (Color on-line.) 

interacting systems (notice how the U = 2 response lies 
above the U = 1.5 and U = 1 responses for long times). 
We will see this phenomenon recur as we increase the 
electric field further. As mentioned above, the accuracy 
continues to improve as we increase the field. For all 
cases shown, the error in the zeroth and second moments 
is less than 1% for times larger than 5. 

In Fig. [5l we plot results for the current with E = 1. 
In panel (a) we show normal metals [/ = 0, 0.125, 0.25, 
and 0.5, while in panel (b) we show one anomalous metal 
U = 1 and three Mott insulators {U = 1.5, 2, and 3). 
Note how the metals continue to show damped Bloch 
oscillations, and that we cannot reach the steady state 
within the time window. It may appear that the U = 0.5 
case is not oscillating with the Bloch period, but that de- 
viation is most likely coming from some higher harmonics 
in the transient response that should vanish as time in- 
creases, because the oscillations are clearly not with some 
slightly larger period in this range. The accuracy for all 
of these metallic cases is better than 1% for both the 
zeroth and second moment for times greater than t = 5. 

Panel (b) shows the current for Mott insulators (and 
one anomalous metal). As the system becomes more 
Mott-insulating, we see more irregular oscillations of a 




Time 

FIG. 5: Scaled nonequilibrium current for different values of 
U with E = 1. In panel (a), we show cases that are metals 
in equilibrium {U — 0, 0.125, 0.25, and 0.5. In panel (b), we 
show one anomalous metal {U = 1) and three Mott insulators 
{U = 1.5, 2 and 3); note the reduced scale for the vertical axis 
versus that in panel (a). We used a quadratic extrapolation 
formula for U < 0.5 {At = 0.1, 0.067, and 0.05), a cubic 
extrapolation formula for U = 1 {At = 0.1, 0.067, 0.05, and 
0.04), a quadratic extrapolation for U — 1.5 {At = 0.067, 
0.05, and 0.04), a quartic extrapolation formula for U = 2 
{At = 0.067, 0.05, 0.04, 0.033, and 0.025), and a quadratic 
extrapolation formula for (7 = 3 {At — 0.02, 0.0167, and 
0.0143). (Color on-line.) 



larger amplitude. Note how, in addition, the C/ = 3 Mott- 
insulating current lies above the U = 2 Mott-insulating 
current at long times. For smaller U values, the oscil- 
lations still maintain a fair amount of regularity with a 
period close to the Bloch period. As shown in Ref. 0, the 
irregular oscillations continue out to long times. These 
results clearly show the emerging trends as we move from 
a metal to a Mott insulator. The current initially is 
quenched via a damping of the Bloch oscillations, but as 
the coupling strength increases, the oscillations become 
more irregular and lower in overall amplitude. As the 
field strength increases, we see the oscillations survive 
out to longer and longer times. In most cases, we are 
unable to reach the steady-state limit within our time 
window. The accuracy of our calculations is still quite 
good. The extrapolated results have errors of less than 
1% for for both moments when U = 1, less than 2% for 
both moments when U = 1.5 and U = 2, and less than 
3% for both moments when U = 3. 

In Fig. [51 we plot the current for our final electric field 
E = 2. The behavior in this case is quite different from 
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FIG. 6: Scaled nonequilibrium current for different values of 
U with E = 2. In panel (a), we show cases that are metals 
in equihbrium {U = 0, 0.125, 0.25, 0.5 and 1). In panel (b), 
we show Mott insulators {U = 1.5, 2, 3, and 4); note the 
reduced scale for the vertical axis versus that in panel (a). 
We could not fit the U labels onto panel (a). The U = 
case (black) has no damping to the oscillations. The U = 
0.125 case (red) looks like a damped oscillator, because the 
beat period is about 50. The U = 0.25 case (green) shows 
one beat period (it has the second largest amplitude near 
t — 30). The U — 0.5 case (blue) shows a few beats and 
is the third largest amplitude near t = 30. The U — 1 case 
(magenta) is the most irregular looking of the curves. We used 
a quadratic extrapolation formula U < 0.5 {At = 0.1, 0.067, 
and 0.05), a cubic extrapolation formula for U = 1 (At = 0.1, 
0.067, 0.05, and 0.04), and for U = 1.5 (At = 0.067, 0.05, 
0.04, and 0.033), a quadratic extrapolation formula for U = 2 
{At — 0.04, 0.033, and 0.025), a cubic extrapolation formula 
for f/ = 3 {At = 0.04, 0.033, 0.025, and 0.0167), and a linear 
extrapolation ior U = i {At = 0.0167 and 0.0143). (Color 
on-line.) 



the other cases. In the metallic case, we do not have 
simple damped oscillations. Instead, the oscillations ini- 
tially decay but then grow, reminiscent of damped beats. 
The beat period turns out to be 2tt/U, decreasing as 
U increases. This result was predicted from nonequilib- 
rium calculations for two particles on a one-dimensional 
chainS^, but we see here that it appears to only enter 
once the electric field is large enough (for an infinite- 
dimensional lattice). We will treat this phenomenon in 
more detail below and resolve the underlying physical 
mechanism behind its occurrence. The beat period be- 
comes so short, that the U = 1 results look like quite 
irregular oscillations. 

The beats disappear as we pass through the metal- 
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FIG. 7: Unsealed nonequilibrium current (At = 0.1) for 
S = 1 and (a) U = 0.125, (b) U = 0.25, and (c) U = 0.5. The 
longer time cutoff allows us to see the initial development of 
beat-like behavior in the current response. 



insulator transition. Instead, we see the irregularity in 
the oscillations become more apparent. In panel (b), one 
can see that the near-critical Mott insulator {U — 1.5) 
still looks like it could be respresented in terms of damped 
oscillations with beats, but the oscillations are not sim- 
ple sinusoidal oscillations. As U is increased further, the 
oscillations are sharply attenuated for times larger than 
15. They continue to have significant short-period irreg- 
ular oscillations, which appear to develop more strongly 
as U increases. In addition, we see all of the different 
Mott insulators to have small current at long times, and 
we do not see some cases having significantly larger cur- 
rent as we did when E — 1. In spite of our ability to 
extrapolate our numerical calculations to At 0, these 
Mott-insulating results are the most difficult to obtain re- 
liable data for, and it is quite likely that the final curves 
for the current have pointwise errors on the order of 10%, 
but we have no way of rigorously estimating the errors, 
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because our calculations have been pushed to the limit of 
what we can reliable achieve. The errors for the metallic 
cases are all quite low (less than 1% for both moments 
for all metals and for U = 1.5). But they grow in the 
Mott insulators (less than 5% error for the first moment, 
and less than 7% error for the second moment for larger 
U values). 
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FIG. 8: Unsealed nonequilibrium current (At = 0.1) for 
E = 2imd (a) U = 0.125, (b) U = 0.25, and (c) U = 0.5. Here 
the beats are clearly visible, but one can also see a dephasing 
at long times, because the beat amplitude does not go to zero. 

Because the phenomena of beats newly develops for 
large electric fields in metals, we need to investigate the 
origin of the beats more thoroughly. To do this, we per- 
form a series of calculations with At — 0.1 and no scaling 
of the data, but with a much longer time window {t runs 
out to 195). We ran calculations for E = 0.5, E = 1, 
and E = 2. In the case when E = 0.5, we saw no indica- 
tion of any beat phenomena in the data. The results for 
E — 1 are shown in Fig. [71 Panel (a) is for U = 0.125, 
panel (b) for U = 0.25, and panel (c) for U = 0.5. Note 
how there is clearly some behavior that is reminiscent 



of beats in this data (especially for small U), but it is 
not well formed beats yet. [Note that the increase in 
amplitude of the current at long times in panel (c) is a 
systematic error associated with the large discretization 
size for the U value and the boundary in time.] We do 
see the period of the beat-like phenomena to decrease as 
U increases, and it appears to be equal to 27r//7. When 
we move to the larger field oi U = 2, which is plotted in 
Fig. [H the beats become completely transparent, even if 
the amplitude modulation becomes somewhat dephased 
at long times. We clearly see the beat period decreasing 
as U increases and being inversely proportional to U, 
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FIG. 9: Local many-body density of states with U — 0.5 
for time T — 95 units after the electric field was turned on 
and (a) E — 0.5, (b) E = 1, and (c) -E = 2; the data are 
from unsealed results with At — 0.1. The Wannier-Strak lad- 
der would have delta functions at the Bloeh frequencies nE. 
Here the many-body density of states is broadened into bands 
about the Bloeh frequencies with a width of approximately 
U. Note how the central band becomes sharply peaked at 
u — ±U /2 as E increases. This is the source of the two fre- 
quencies separated by U that creates the beats in the current 
as a function of time. 
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One can ask what the cause of these beats can be? 
Since the beat period is 27r/[/, we immediately suspect 
the result is coming from two frequencies equal to the 
Bloch frequency plus or minus J7/2. Then the combina- 
tion of these oscillations will produce the expected beats. 
If we extract the retarded Green's function from the 
contour-ordered Green's function, transform to Wigner 
coordinates and perform a Fourier transform with respect 
to the relative time, we can then plot the local many- 
body density of states as a function of average time. For 
large times, the density of states approaches a steady 
state limit, and we plot a large time result (T = 95) in 
Fig. [9] for U = 0.5 and the three different E values (0.5, 
1, and 2). One can see that as the electric field is in- 
creased, the DOS develops two sharp peaks centered at 
uj = zLU/2. While we can see these peaks already for 
E = 1, they really become well separated and distinct 
for i? — 2. We believe that these two delta-function- 
like peaks, when combined with the underlying Bloch 
frequency, are the source of the beat phenomena in the 
current, and since we do not see the double peak struc- 
ture develop for smaller fields, the beat behavior in the 
current is a strong field behavior. This is hinted at in 
the work with two particles in one-dimensionS^, but that 
work showed the additional oscillations occurring for all 
U and E. As one increases the particle density, and the 
spatial dimension, the occurence of beats requires a crit- 
ical value of the electric field, which is around E = 1 for 
our case. 



IV. CONCLUSIONS 

In this work, we have shown a detailed derivation of 
the generalization of DMFT to nonequilibrium problems 
given by the so-called Keldysh boundary condition — the 
system starts in equilibrium and then a field is turned 
on and the system is monitored as it is driven to a 
steady state. Our focus was on the current and how the 
Bloch oscillations are quenched and change character as 
electron-electron correlations are increased. We find for 
weak fields the picture in metals is pretty much what one 
might have guessed — the Bloch oscillations are damped 
and the system approaches a steady state with what ap- 
pears to be a constant current (although we cannot rule 
out low amplitude oscillations with the Bloch period su- 
perimposed on top of the constant response). As the field 
is increased, we find the oscillations survive out to longer 
and longer times. When the coupling increases to that 
of a Mott insulator, the oscillations change their char- 
acter and become quite irregular. This process is a slow 
evolution and not a sharp "phase transition" . Most inter- 
esting, is the appearance of beats in the current for large 
electric fields in metals. These beats have a period in- 
versely proportional to the interaction strength and are 



quite robust. As time increases however, the beats do 
dephase, and eventually we get a constant amplitude re- 
sponse (which may be decreasing due to damping) that 
no longer has an oscillatory envelope. We also explained 
in detail many of the technical issues required by the 
discretization and by the numerics to obtain an efficient 
iterative algorithm that scales to many hundreds if not 
thousands of processors. 

While this work is a significant advance in theoretical 
many-body physics, it is not clear that this phenomena 
can be observed directly in solid-state systems. The prob- 
lem there is that most solid state systems where Bloch 
oscillations can be observed have higher energy bands 
present. If the field becomes large enough to induce tun- 
neling between the bands, then the theory we have de- 
rived will not apply, because we have neglected Zener 
tunneling. This implies the observation of beats may 
be difficult to see in a solid-state system. Furthermore, 
the fields are so large here that the Bloch frequencies 
are enormous, with periods probably lying in the fem- 
tosecond range. It is a challenge to find experimental 
techniques to measure such rapid oscillations. 

One system that might be appropriate for this exper- 
iment, however, is ultracold atoms placed in optical lat- 
tices. If we consider a mixture of two species where one of 
the atomic species is much more massive than the other, 
and if the delocalized fermions are in a band that is well 
separated from the higher-energy bands, then by detun- 
ing the optical lattice, so that it appears as a static lattice 
in a moving frame (via the Doppler effect), one can "pull" 
the lattice through the atomic cloud, in direct analogy to 
the Hamiltonian gauge calculations presented here. Un- 
fortunately there is no direct way measure the current, 
but it could be reconstructed, in principle, from a time- 
of-flight measurement which determined the distribution 
of atoms through the Brillouin zone. 

In future work, we plan to examine the steady-state 
limit directly and how the distribution functions tran- 
siently evolve after the field is turned on (the latter being 
important for direct comparison to time-of-flight mea- 
surements in ultracold atoms). 
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